18  Production model - river 4 (IS)

Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.

In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv

Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R

There are 3 models here: All cohorts:
model 1, with quadratic weight effect on phi By cohort:
model 2, with linear weight effect on phi model 3, with quadratic weight effect on phi

18.0.1 How many ageInSamples to include?

Code
all <- tar_read(cdWB_electro_target)
table(all %>% filter(river == "wb obear")|> dplyr::select(ageInSamples))
ageInSamples
   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
1308 1890 1007  837 1486  875  446  298  444  159  155   57   47   17   16    3 
  16   18 
   3    3 
Code
cohorts <- 2002:2014

18.1 All Cohorts, modelNum 1: phi(i,t,c) * g(i,t,c), p(i,t) model

18.2 By Cohort, modelNum 2: phi(i,t) * g(i,t), p(i,t) model

Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

      logit(phi[i,t]) <- 
        phiInt[i,t] + 
        phiBeta[1,i,t] * weight[i,t] 

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t)

Survival/growth rate interaction:

Multiplicative

Model runs:

18.2.1 Retrieve model results

Code
library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
# 

# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
source('./models/production/modelProduction_OB_functionsToSource.R')
source('./models/production/qmdProduction_OB_functionsToSource.R')

modelNum <- 2 # phi * growth

18.2.2 Model code

Code
# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_2_2002_mostRecent.RData')
d$modelCode
{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:1) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in 1:T) {
            weightZ01[i, t] <- weightIfAvailable(weight[i, t], 
                sdWeight, meanWeight, z[i, t], available[i, t])
            z01[i, t] <- aliveIfAvailable(z[i, t], available[i, 
                t])
        }
    }
    for (t in 1:T) {
        sumWeightZ01[t] <- sum(weightZ01[1:N, t])
        sumZ01[t] <- sum(z01[1:N, t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}
Code
getSummaryByCohort <- function(cohort, modelNum = 3) {
  # Load the data
  load(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '_mostRecent.RData')) 
  
  # Input data
  eh <- tar_read_raw(paste0('eh_OB_', cohort, '_target'))
  
  print(paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")))

  MCMCplot(object = d$mcmc, params = c("sumZ01"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("sumWeightZ01"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE)
  MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE)

  get_MCMCTrace_OB(d)
  
   # Create the summary list
  summary_list <- list(
    eh = eh,
    stats = kable(as_tibble(d$runData), caption = "Run statistics"),
    runTime = paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")),
    #traces = get_MCMCTrace_OB(d),
    summaryMod3_tt_growth = MCMCsummary(object = d$mcmc, params = c(
        #"phiIntT", 
        "phiIntTOut", 
        "p", 
        "grIntT", 
        "phiBetaT",
        "sumWeightZ01",
        "sumZ01"
      ), round = 3) %>%
      rownames_to_column(var = "var"),
    summary_OB = get_summary_OB(d, eh)
  )
  
  return(summary_list)
}

#s=summary_OB[, order(summary_OB$ind, summary_OB$t)]

18.2.3 Summaries by cohort

Code
s2002_2 = getSummaryByCohort(2002, modelNum)
[1] "Run time = 21.474 mins"

Code
s2003_2 = getSummaryByCohort(2003, modelNum)
[1] "Run time = 27.403 mins"

Code
s2004_2 = getSummaryByCohort(2004, modelNum)
[1] "Run time = 17.544 mins"

Code
s2005_2 = getSummaryByCohort(2005, modelNum)
[1] "Run time = 8.699 mins"

Code
s2006_2 = getSummaryByCohort(2006, modelNum)
[1] "Run time = 14.192 mins"

Code
s2007_2 = getSummaryByCohort(2007, modelNum)
[1] "Run time = 11.018 mins"

Code
s2008_2 = getSummaryByCohort(2008, modelNum)
[1] "Run time = 9.901 mins"

Code
s2009_2 = getSummaryByCohort(2009, modelNum)
[1] "Run time = 55.498 mins"

Code
s2010_2 = getSummaryByCohort(2010, modelNum)
[1] "Run time = 7.013 mins"

Code
s2011_2 = getSummaryByCohort(2011, modelNum)
[1] "Run time = 6.584 mins"

Code
s2012_2 = getSummaryByCohort(2012, modelNum)
[1] "Run time = 1.09 hours"

Code
s2013_2 = getSummaryByCohort(2013, modelNum)
[1] "Run time = 8.649 mins"

Code
s2014_2 = getSummaryByCohort(2014, modelNum)
[1] "Run time = 6.429 mins"

18.2.4 Combine cohort summaries

Code
sAll_2 <- tibble(.rows = 0) |>
  add_column(!!!s2002_2$summary_OB[0,])

# Fill tibble
for (cohort in cohorts) {
  summary_obj_2 <- get(paste0("s", cohort, "_2"))$summary_OB
  summary_obj_2$cohort <- cohort
  sAll_2 <- bind_rows(sAll_2, summary_obj_2)
}

sAll_2$cohort_ind <- paste0(sAll_2$cohort, "_", sAll_2$ind)
Code
#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_2 = sAll_2)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

18.2.4.1 Select cohort(s):

Code
summary_OB_OJS_all_T_2 = transpose(summary_OB_OJS_all_2)
Code
import { range } from "d3";

viewof selected_cohort_2 = Inputs.select(d3.range(2002, 2015), {
  label: "Which cohort?",
  value: 2002,
  step: 1,
  multiple: true
})
Code
summary_OB_OJS_all_T_selected_cohort_2 = summary_OB_OJS_all_T_2.filter((d) =>
  (selected_cohort_2.includes(d.cohort))
)

18.2.4.2 Select one or more individuals:

Code
uniqueInds_2 = [...new Set(summary_OB_OJS_all_T_selected_cohort_2.map(d => d.cohort_ind))].sort();

viewof selectInd_mod_2 = Inputs.select(uniqueInds_2, {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})


summary_OB_OJS_all_T_selected_cohort_ind_2 = summary_OB_OJS_all_T_selected_cohort_2.filter((d) =>
  (selectInd_mod_2.includes(d.cohort_ind))
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

18.2.4.3 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Probability of survival" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "pSurv"
      }),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "enc8",
        fill: "orange"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "pSurv"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "first",
        y: 1,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "last",
        y: 1,
        "stroke": "red"
      })
      ,
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2, 
                Plot.selectFirst({
                  x: 11,
                  y: 1,
                  text: d => d.cohort
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_2,
      x: "ind"
    }
  })
Code
Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Standardized body mass (mg)" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "mean_weight"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "t",
        y: "mean_weight"
      }),
      <!-- Plot.dot(d, { -->
          <!--   x: "t", -->
          <!--   y: "mean_gr", -->
          <!--   fill: "green" -->
          <!-- }), -->
        Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
          x: "t",
          y: "weightDATA_std",
          fill: "orange"
        }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "first",
        y: 3,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
        x: "last",
        y: 3,
        "stroke": "red"
      }),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
                Plot.selectFirst({
                  x: 1,
                  y: 4,
                  frameAnchor: "top-left",
                  text: (d) => "Cohort = " + d.cohort
                })
      ),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
                Plot.selectFirst({
                  x: 1,
                  y: 3.5,
                  frameAnchor: "top-left",
                  text: (d) => "Residual = " + d.meanResid
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_2,
      x: "ind"
    }
  })

18.2.5 Compare cohorts

Combine scohort$summaryMod2_tt_growth data into one data frame

Code
sParamsAll0_2 <- tibble(.rows = 0) |>
  add_column(!!!s2002_2$summaryMod3_tt_growth[0,]) # need to change to mod2

# Fill tibble
for (cohort in cohorts) {
  summary_obj_2 <- get(paste0("s", cohort, "_2"))$summaryMod3_tt_growth # need to change to mod2
  summary_obj_2$cohort <- cohort
  sParamsAll0_2 <- bind_rows(sParamsAll0_2, summary_obj_2)
}

# remove brackets for filtering
sParamsAll0_2$varNoIndex <- sParamsAll0_2$var |>
  str_remove("\\[.*\\]")

# for numeric ordering
sParamsAll1_2 <- sParamsAll0_2 |> filter(varNoIndex != "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+)\\]") %>%
      str_remove_all("\\[|\\]") %>%
      as.numeric()
  )

sParamsAll2_2 <- sParamsAll0_2 |> filter(varNoIndex == "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    varNumeric2 = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric()
  )

sParamsAll_2 = bind_rows(sParamsAll1_2, sParamsAll2_2)

18.2.5.1 Combined cohorts

18.2.5.2 p

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.5.3 phi

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.5.4 grIntT

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.5.5 grIntT, filtered

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.5.6 phiBetaT

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples") +
  facet_wrap((~varNumeric1))

18.2.5.7 sumWeightZ01

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.5.8 sumZ01

Code
ggplot(sParamsAll_2 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.2.6 Size-dependent survival

Need to get rid of squared term for cohort-specific models

Code
phiInt_2 <- sParamsAll_2 |> filter(varNoIndex == "phiIntTOut") |>
  dplyr::select(mean, varNumeric1, cohort) |>
  rename(phiIntOut = mean, ais = varNumeric1)

phiBetas_2 <- sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
  dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
  pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
  rename(ais = varNumeric2)

phiSize_2 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |> 
   left_join(phiInt_2) |> 
   left_join(phiBetas_2) |>
   mutate(pSurv = ilogit(phiIntOut + beta1*weight))# + beta_2*weight^2))

ggplot(phiSize_2, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ais)

18.3 By Cohort, modelNum 3: phi(i,t) * g(i,t)^2, p(i,t) model

Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.

Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)

Probability of survival (phi) model structure:

      logit(phi[i,t]) <- 
        phiInt[i,t] + 
        phiBeta[1,i,t] * weight[i,t] + 
        phiBeta[2,i,t] * weight[i,t] * weight[i,t] 

Probability of capture (p) model structure:

p(t,i) <- pInt(t,i)

Growth rate (gr) model structure:

gr(t,i) <- grInt(t)

Survival/growth rate interaction:

Multiplicative

Model runs:

18.3.1 Retrieve model results

Code
library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
# 
# 

# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
#source('./models/production/modelProduction_OB_functionsToSource.R')
#source('./models/production/qmdProduction_OB_functionsToSource.R')

modelNum <- 3 # phi * growth

18.3.2 Model code

Code
# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_3_2002_mostRecent.RData')
d$modelCode
{
    for (i in 1:N) {
        weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]], 
            sd = 2)
        weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]], 
            sd = weightSD)
        for (t in (first[i]):(last[i] - 1)) {
            weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i, 
                t]/90
            weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1], 
                sd = weightSD)
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
        }
    }
    for (t in 1:(T - 1)) {
        grIntT[t] ~ dnorm(0, sd = 1000)
    }
    delta[1] <- 1
    delta[2] <- 0
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            gamma[1, 1, t, i] <- phi[i, t]
            gamma[1, 2, t, i] <- 1 - phi[i, t]
            gamma[2, 1, t, i] <- 0
            gamma[2, 2, t, i] <- 1
        }
    }
    for (t in 1:(T - 1)) {
        p[t] ~ dunif(0, 1)
        omega[1, 1, t] <- 1 - p[t]
        omega[1, 2, t] <- p[t]
        omega[2, 1, t] <- 1
        omega[2, 2, t] <- 0
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i, 
                t] * weight[i, t] + phiBeta[2, i, t] * weight[i, 
                t] * weight[i, t]
        }
    }
    for (i in 1:N) {
        for (t in first[i]:(last[i] - 1)) {
            phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
            for (b in 1:2) {
                phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
            }
        }
    }
    for (t in 1:(T - 1)) {
        phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
        phiIntTOut[t] <- ilogit(phiIntT[t])
        phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
        phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
    }
    for (i in 1:N) {
        for (t in 1:T) {
            weightZ01[i, t] <- weightIfAvailable(weight[i, t], 
                sdWeight, meanWeight, z[i, t], available[i, t])
            z01[i, t] <- aliveIfAvailable(z[i, t], available[i, 
                t])
        }
    }
    for (t in 1:T) {
        sumWeightZ01[t] <- sum(weightZ01[1:N, t])
        sumZ01[t] <- sum(z01[1:N, t])
    }
    for (i in 1:N) {
        z[i, first[i]] ~ dcat(delta[1:2])
        for (j in first[i]:(last[i] - 1)) {
            z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
            yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
        }
    }
}

18.3.3 Summaries by cohort

Code
s2002_3 = getSummaryByCohort(2002, modelNum)
[1] "Run time = 32.092 mins"

Code
s2003_3 = getSummaryByCohort(2003, modelNum)
[1] "Run time = 37.935 mins"

Code
s2004_3 = getSummaryByCohort(2004, modelNum)
[1] "Run time = 28.262 mins"

Code
s2005_3 = getSummaryByCohort(2005, modelNum)
[1] "Run time = 12.692 mins"

Code
s2006_3 = getSummaryByCohort(2006, modelNum)
[1] "Run time = 20.456 mins"

Code
s2007_3 = getSummaryByCohort(2007, modelNum)
[1] "Run time = 15.911 mins"

Code
s2008_3 = getSummaryByCohort(2008, modelNum)
[1] "Run time = 14.706 mins"

Code
s2009_3 = getSummaryByCohort(2009, modelNum)
[1] "Run time = 1.453 hours"

Code
s2010_3 = getSummaryByCohort(2010, modelNum)
[1] "Run time = 10.1 mins"

Code
s2011_3 = getSummaryByCohort(2011, modelNum)
[1] "Run time = 9.941 mins"

Code
s2012_3 = getSummaryByCohort(2012, modelNum)
[1] "Run time = 1.39 hours"

Code
s2013_3 = getSummaryByCohort(2013, modelNum)
[1] "Run time = 13.827 mins"

Code
s2014_3 = getSummaryByCohort(2014, modelNum)
[1] "Run time = 8.986 mins"

18.3.4 Combine cohort summaries

Code
sAll_3 <- tibble(.rows = 0) |>
  add_column(!!!s2002_3$summary_OB[0,])

# Fill tibble
for (cohort in cohorts) {
  summary_obj_3 <- get(paste0("s", cohort, "_3"))$summary_OB
  summary_obj_3$cohort <- cohort
  sAll_3 <- bind_rows(sAll_3, summary_obj_3)
}

sAll_3$cohort_ind <- paste0(sAll_3$cohort, "_", sAll_3$ind)
Code
#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_3 = sAll_3)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))

18.3.4.1 Select cohort(s):

Code
summary_OB_OJS_all_T_3 = transpose(summary_OB_OJS_all_3)
Code
import { range } from "d3";

viewof selected_cohort_3 = Inputs.select(d3.range(2002, 2015), {
  label: "Which cohort?",
  value: 2002,
  step: 1,
  multiple: true
})
Code
summary_OB_OJS_all_T_selected_cohort_3 = summary_OB_OJS_all_T_3.filter((d) =>
  (selected_cohort_3.includes(d.cohort))
)

18.3.4.2 Select one or more individuals:

Code
uniqueInds_3 = [...new Set(summary_OB_OJS_all_T_selected_cohort_3.map(d => d.cohort_ind))].sort();

viewof selectInd_mod_3 = Inputs.select(uniqueInds_3, {
  label: "Which fish?",
  value: 1,
  step: 1,
  multiple: true
})


summary_OB_OJS_all_T_selected_cohort_ind_3 = summary_OB_OJS_all_T_selected_cohort_3.filter((d) =>
  (selectInd_mod_3.includes(d.cohort_ind))
)

Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

18.3.4.3 Plot survival

Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.

Code
Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Probability of survival" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "pSurv"
      }),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "enc8",
        fill: "orange"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "pSurv"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "first",
        y: 1,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "last",
        y: 1,
        "stroke": "red"
      })
      ,
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3, 
                Plot.selectFirst({
                  x: 11,
                  y: 1,
                  text: d => d.cohort
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_3,
      x: "ind"
    }
  })
Code
Plot.plot({
    width: width,
    height: 350,
    inset: 10,
    color: {
      scheme: "greys"
    },
    x: { label: "Age/season combination" },
    y: { label: "Standardized body mass (mg)" },
    marks: [
      Plot.frame(),
      Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "mean_weight"
      }),
      Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "t",
        y: "mean_weight"
      }),
      <!-- Plot.dot(d, { -->
          <!--   x: "t", -->
          <!--   y: "mean_gr", -->
          <!--   fill: "green" -->
          <!-- }), -->
        Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
          x: "t",
          y: "weightDATA_std",
          fill: "orange"
        }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "first",
        y: 3,
        "stroke": "green"
      }),
      Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
        x: "last",
        y: 3,
        "stroke": "red"
      }),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
                Plot.selectFirst({
                  x: 1,
                  y: 4,
                  frameAnchor: "top-left",
                  text: (d) => "Cohort = " + d.cohort
                })
      ),
      Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
                Plot.selectFirst({
                  x: 1,
                  y: 3.5,
                  frameAnchor: "top-left",
                  text: (d) => "Residual = " + d.meanResid
                })
      )
    ],
    facet: {
      data: summary_OB_OJS_all_T_selected_cohort_ind_3,
      x: "ind"
    }
  })

18.3.5 Compare cohorts

Combine scohort$summaryMod2_tt_growth data into one data frame

Code
sParamsAll0_3 <- tibble(.rows = 0) |>
  add_column(!!!s2002_3$summaryMod3_tt_growth[0,]) 

# Fill tibble
for (cohort in cohorts) {
  summary_obj_3 <- get(paste0("s", cohort, "_3"))$summaryMod3_tt_growth 
  summary_obj_3$cohort <- cohort
  sParamsAll0_3 <- bind_rows(sParamsAll0_3, summary_obj_3)
}

# remove brackets for filtering
sParamsAll0_3$varNoIndex <- sParamsAll0_3$var |>
  str_remove("\\[.*\\]")

# for numeric ordering
sParamsAll1_3 <- sParamsAll0_3 |> filter(varNoIndex != "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+)\\]") %>%
      str_remove_all("\\[|\\]") %>%
      as.numeric()
  )

sParamsAll2_3 <- sParamsAll0_3 |> filter(varNoIndex == "phiBetaT") |>
  mutate(
    varNumeric1 = var |>
      str_extract("\\[(\\d+),") %>%           # Changed to match number before comma
      str_remove_all("\\[|,") %>%             # Remove [ and comma
      as.numeric(),
    varNumeric2 = var |>
      str_extract(",\\s*(\\d+)\\]") %>%
      str_remove_all(",|\\s|\\]") %>%
      as.numeric()
  )

sParamsAll_3 = bind_rows(sParamsAll1_3, sParamsAll2_3)

18.3.5.1 Combined cohorts

18.3.5.2 p

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.5.3 phi

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.5.4 grIntT

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.5.5 grIntT, filtered

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.5.6 phiBetaT

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples") +
  facet_wrap((~varNumeric1))

18.3.5.7 sumWeightZ01

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.5.8 sumZ01

Code
ggplot(sParamsAll_3 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
  geom_point() + 
  geom_line() +
  xlab("AgeInSamples")

18.3.6 Size-dependent survival

Need to get rid of squared term for cohort-specific models

Code
phiInt_3 <- sParamsAll_3 |> filter(varNoIndex == "phiIntTOut") |>
  dplyr::select(mean, varNumeric1, cohort) |>
  rename(phiIntOut = mean, ais = varNumeric1)

phiBetas_3 <- sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
  dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
  pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
  rename(ais = varNumeric2)

phiSize_3 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |> 
   left_join(phiInt_3) |> 
   left_join(phiBetas_3) |>
   mutate(pSurv = ilogit(phiIntOut + beta1*weight + beta2*weight^2))

ggplot(phiSize_3, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ais)